In zozlak/WicPomiarSatelitarny: Pakiet ulatwiajacy przeprowadzenie warsztatow z pomiaru satelitarnego

Przygotowanie

Załadujmy potrzebne pakiety.

library(WicPomiarSatelitarny)

Puławy

Pobranie danych

Pobierzmy zdjęcia od 29 czerwca do dziś, ale tylko takie, dla których zachmurzenie nie przekraczało 60%.
(pozycja Puław to z grubsza 21.95 E, 51.5 N)

zdjPulawy = rbind(
  pobierzZdjecia(21.95, 51.5, "2016-06-29", "2016-06-29", 60, "dane/Pulawy", "wget"),
  pobierzZdjecia(21.95, 51.5, "2016-08-08", "2016-08-08", 60, "dane/Pulawy", "wget")
)

Generowanie obrazów kolorowych

# stwórz listę obiektów R na podstawie zdjęć zapisanych na dysku
zdjecia = wczytajZdjecia("dane/Pulawy", "2016-08-08")

# skomponuj obraz kolorowy z trzech wskazanych pasm
eksportujKolor(zdjecia[["B04"]], zdjecia[["B03"]], zdjecia[["B02"]], 'dane/Pulawy/rgb.tif', 0, 5000)

# pokoloruj klasyfikację scen w zadany sposób
eksportujMono(zdjecia[["SCL"]], "dane/Pulawy/scl.tif", 0, 11, c('red', 'red', 'white', 'white', 'green', 'brown', 'blue', 'white', 'white', 'white', 'white', 'white'))

# pokoloruj pasmo wegetacji na zielono
eksportujMono(zdjecia[["B8A"]], "dane/Pulawy/B8A_kolor.tif", 0, 5000, c('black', 'green'))

NDVI

Policzmy wskaźnik NDVI.

zdjecia[["NDVI"]] = (zdjecia[["B8A"]] - zdjecia[["B04"]]) / (zdjecia[["B8A"]] + zdjecia[["B04"]])
eksportujMono(zdjecia[["NDVI"]], "dane/Pulawy/NDVI_2016-08-08.tif", -1, 1, c('blue', 'black', 'green'))

Analogicznie policzmy NDVI dla 2016-06-29

zdjecia2 = wczytajZdjecia("dane/Pulawy", "2016-06-29")
zdjecia2[["NDVI"]] = (zdjecia2[["B8A"]] - zdjecia2[["B04"]]) / (zdjecia2[["B8A"]] + zdjecia2[["B04"]])
eksportujMono(zdjecia2[["NDVI"]], "dane/Pulawy/NDVI_2016-06-29.tif", -1, 1, c('blue', 'black', 'green'))

Chmury, itp.

Ustawmy na braki danych wszystkie piksele niesklasyfikowane jako "roślinność", "goła ziemia" ani "woda" (liczenie NDVI dla chmur, braków danych, itp. nie ma sensu).

zdjecia[["NDVI"]][! zdjecia[["SCL"]] == 4] = NA
zdjecia2[["NDVI"]][! zdjecia2[["SCL"]] %in% c(4, 5, 6)] = NA
eksportujMono(zdjecia[["NDVI"]], "dane/Pulawy/NDVI_bez_chmur_2016-08-08.tif", -1, 1, c('blue', 'black', 'green'))
eksportujMono(zdjecia2[["NDVI"]], "dane/Pulawy/NDVI_bez_chmur_2016-06-29.tif", -1, 1, c('blue', 'black', 'green'))

NDVI po raz drugi

Oznaczmy miejsca, wg różnicy wartości NDVI pomiędzy 29 czerwca i 8 sierpnia.

roznica = zdjecia[["NDVI"]] - zdjecia2[["NDVI"]]
eksportujMono(roznica, "dane/Pulawy/roznica.tif", -1, 1, c('red', 'black', 'green'))

McMurray

W maju 2016 r. w okolicach miasta McMurray w Kanadzie miał miejsce olbrzymi pożar lasu (trwał ponad miesiąc). Porównując NDVI sprzed i po pożarze spróbujemy odnaleźć pogorzeliska.

Pobieranie danych

Współrzędne miasteczka to z grubsza 111.37 W, 56.7 N. Pobierzemy dane z 12 kwietnia i 21 czerwca.

zdjMcMurray = rbind(
  pobierzZdjecia(-111.37, 56.7, "2016-04-12", "2016-04-12", 60, "dane/Pulawy", "wget"),
  pobierzZdjecia(-111.37, 56.7, "2016-06-21", "2016-06-21", 60, "dane/Pulawy", "wget")
)

Generowanie obrazów

Wczytajmy obrazy.

przed = wczytajZdjecia("dane/McMurray/", "2016-04-12")
po = wczytajZdjecia("dane/McMurray/", "2016-06-21")

Policzmy NDVI i jego różnicę.

ndviPrzed = (przed[["B8A"]] - przed[["B04"]]) / (przed[["B8A"]] + przed[["B04"]])
ndviPo = (po[["B8A"]] - po[["B04"]]) / (po[["B8A"]] + po[["B04"]])
roznica = ndviPo - ndviPrzed

Za pomocą klasyfikacji scen zamieńmy na braki danych piksele przysłonięte chmurami na dowolnym z obrazów (różnica wartości NDVI dla tych pikseli nic nam nie powie).

roznica2[! przed[["SCL"]] %in% 4:6 | ! po[["SCL"]] %in% 4:6] = NA

Wyeksportujmy obrazki.

Obszary oznaczone na czerwono na obrazie "roznica2.tiff" w północnej części zdjęcia to prawdopodobnie pogorzeliska.

eksportujMono(roznica, 'dane/McMurray/roznica.tiff', -1, 1, c('red', 'black', 'green'))
eksportujMono(roznica2, 'dane/McMurray/roznica2.tiff', -1, 1, c('red', 'black', 'green'))

eksportujMono(ndviPrzed, 'dane/McMurray/ndvi_przed.tiff', -1, 1, c('blue', 'black', 'green'))
eksportujMono(ndviPo, 'dane/McMurray/ndvi_po.tiff', -1, 1, c('blue', 'black', 'green'))

eksportujMono(przed[["SCL"]], 'dane/McMurray/scl_przed.tiff', 0, 11, c('red', 'red', 'white', 'white', 'green', 'brown', 'blue', 'white', 'white', 'white', 'white', 'white'))
eksportujMono(po[["SCL"]], 'dane/McMurray/scl_po.tiff', 0, 11, c('red', 'red', 'white', 'white', 'green', 'brown', 'blue', 'white', 'white', 'white', 'white', 'white'))

eksportujKolor(przed[["B04"]], przed[["B03"]], przed[["B02"]], 'dane/McMurray/przed_rgb.tiff', 0, 5000)
eksportujKolor(po[["B04"]], po[["B03"]], po[["B02"]], 'dane/McMurray/po_rgb.tiff', 0, 5000)